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I.  INTRODUCTION 


Estimating  the  power  spectrum  associated  with  a  2-D  random  process  is  important 
in  a  number  of  applications.  In  digital  image  processing,  for  example,  Wiener  filtering 
may  be  used  to  solve  image  restoration  problems  in  which  a  signal  is  degraded  by 
additive  random  noise  [Ref.  1].  The  frequency  response  of  a  noncausal  Wiener 
filter  requires  knowledge  of  the  spectral  contents  of  the  signal  and  background  noise. 
Another  example  is  an  array  of  sensors.  The  jfrequency  wave  number  spectrum  generated 
from  such  an  array  contains  information  about  the  signal  source  and  direction  of  arrival. 

Various  methods  of  spectral  estimation  are  discussed  in  the  literature  [Refs. 
1,  2,  3].  This  tf^esis  treats  the  topic  of  autoregressive  (AR)  spectral  estimation. 
Because  this  method  is  known  to  produce  high  resolution  spectral  estimates  with  a  small 
data  set,  it  has  attracted  considerable  interest  [Refs.  4,  5,  6,  7]. 

Determining  estimates  of  AR  parameters  requires  the  solution  of  a  set  of  normal 
equations.  This  in  turn  requires  the  inversion  of  an  AxTV  autocorrelation  matrix.  Direct 
inversion  of  the  autocorrelation  matrix  requires  0{N^)  multiplications,  which  is 
computationally  intensive  for  large  N.  An  iterative  method  is  presented  in  this  thesis 
which  reduces  the  number  of  multiplications  significantly. 

This  thesis  is  organized  into  five  chapters.  The  remainder  of  this  chapter  is  devoted 
to  defining  the  power  spectral  density  and  the  resulting  forms  of  the  normal  equation 
when  quarter-plane  (QP)  and  nonsymmetric  half-plane  (NSHP)  support  are  used.  Chapter 
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II  develops  an  iterative  solution  to  the  normal  equation  for  both  QP  and  NSHP  support 
which  assumes  the  autocorrelation  matrix  is  toeplitz  block-tocplitz.  Chapter  ni  discusses 
the  use  of  the  covariance  method  to  estimate  the  autocorrelation  matrix.  Because  the 
resulting  autocorrelation  matrix  for  the  covariance  method  is  not  toeplitz  block-tocplitz, 
a  modification  of  the  iterative  solution  presented  in  Chapter  11  is  developed. 
Additionally,  the  technique  of  combined-quadrant  (CQ)  averaging  is  introduced  in  this 
chapter.  Results  obtained  from  using  QP  and  NSHP  support,  as  well  as  CQ  averaging,  to 
estimate  spectral  densities  are  presented  in  Chapter  FV.  Conclusions  and 
recommendations  for  future  study  are  given  in  Chapter  V. 

A.  POWER  SPECTRAL  DENSITY 

A  2-D  random  process  is  a  discrete  function  such  that,  for  each  coordinate 

pair  (//|,«2),  the  value  of  is  a  random  variable.  The  power  spectral  density 

P^(w,,(Oj)  of  a  wide-sense  stationary  random  process  is  the  Fourier  transform  of  the 
autocorrelation  function. 


<10  eo 


(1-1) 


where 


(1-2) 


and  £[•]  denotes  the  statistical  expectation.  Conventional  methods  of  estimating  the 
spectral  density  include  the  correlogram.  In  this  method,  an  estimate  of  the 
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autocorrelation  function  is  substituted  into  (1-1),  where  is  an  estimate 


such  as 


-  -rrrr  E  E  x(n,,n^)x\n^-l^,n^-l^), 


(1-3) 


/=o  /,*o 


and  Nx^2  is  the  number  of  points  in  the  random  process.  The  modem  methods  of  spectral 
estimation  include  representing  the  random  process  with  an  autoregressive  model.  The 
model  parameters  are  then  used  to  estimate  the  spectral  density.  This  is  the  method 
employed  in  this  thesis. 

B.  AUTOREGRESSIVE  SPECTRAL  ESTIMATION 

Autoregressive  spectral  estimation  assumes  the  random  process  Afnpttj)  is  the 

response  of  an  AR  model  excited  by  white  noise  H'fnp/Jj)  with  variance  as  shown  in 

Figure  1.  To  estimate  the  AR  parameters  a(«,,«2),  the  system  is  expressed  as  a  recursive 
difference  equation  given  by 

A(/J,,/J2)  =  -  53  E  +  M'(/7,,/J2)  (1-4) 

(  i.i  )  e  A 


where  A  is  the  region  of  support  over  which  a(i,j)  have  non-zero  values.  The  frequency- 
response  function  of  the  AR  model  is  expressed  as 


//((0p(0  )  -  - - — __ 

eA 


(1-5) 


3 


w(ni,n2) 


1 


x(n  1,(12) 


The  AR  power  spectral  estimate  r^(,(0,,(0j)  is  given  by 


(1-6) 


wliere  is  the  spectral  density  of  the  white  noise  input.  Tlie  input  has  a  constant 


power  spectrum  of  amplitude  0^.  Therefore  (1-6)  can  be  written  as 


r^((o,.(o,)  = 


(1-7) 


To  estimate  the  AR  parajneters,  multiply  both  sides  of  (1-4)  by  x’{n^-l and 
take  the  statistical  expectation  [Ref.  2].  Tliis  leads  to  the  following  nonnal  equation 


4 


(1-8) 


The  structure  of  the  normal  equation  depends  upon  the  shape  of  A.  Two  regions  of 
support  will  be  considered  in  the  following,  quarter-plane  and  nonsymmetric  half-plane 
support. 

1.  Quarter-Plane  Support 

The  region  A  is  said  to  have  quarter-plane  support  when  a{ij)  has  non-zero 
values  in  one  quadrant  only,  as  shown  in  Figure  2,  which  illustrates  quarter  plane  support 
in  the  first  quadrant.  In  this  case  the  normal  equation  becomes 


p,-i  />,-i 

i:  r 

(t.y  )  *  (0.0) 


a{i 


RSI,, I,) 


(1-9) 


where  /,=0,1,2,...,P,-1,  /2=0,1,2,...,F2-1,  and  P,  and  Pj  ^re  the  dimensions  of  /\  as  shown 
in  Figure  2.  If  it  is  assumed  that  a(0,0)  =  1,  then  (1.9)  can  be  written  in  block-matrix 
form  as 


R,  /?, 


" 

" 

- 

R-2  • 

5(0. 

/?_, 

0 

^0 

0 

^P,-3 

0 

(1-10) 
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Figure  2.  QP  stippoil  in  first  quadrant 

wliere  each  block  is  given  by 

R{k,0)  Rik.  \)  • 

^  R{k,l)  R{k.O)  -  R(k,-F,^2) 

R{k,r,-l)  R(k,P^~2)  ••  R{k.d) 

=  [a{k,Q)  a{k,l)  a{k,2)  -  a{k,P^-\)y, 


Note  that  the  matrix  R  has  a  toeplitz  block-toeplitz  structure.  That  is,  the  blocks  along 
the  diagonals  are  equal,  and  the  elements  along  the  diagonals  within  the  blocks  are  also 
equal. 

2.  Nonsymmetric  Half-Plane  Support 

When  a{ij)  has  non-zero  values  in  a  region  of  the  form  shown  in  Figure  3,  the 
region  A  is  referred  to  as  nonsymmetric  half-plane  (NSHP)  sujyport.  In  this  case  the 
normal  equation  becomes 

P,-1  P,-l 

Y.  a (0 ,;)/?,(/,, /j -7)  +  Y.  (1-14) 

/■=l  I'^l  y  = 

(i.J)  5t  (0.0) 


where 

/,  =  0,1,,2,...,P,-1, 

0,1,,2,,...,I,+P,-1,  for/,  =0, 

for  /,  ^  0. 

and  P,  and  Pj  ih®  dimensions  of  A,  and  Lj  Is  a  negative  number  as  defined  in  Figure 
3.  Let  a(0,0)  =  1.  Then  (1-14)  may  be  written  in  matrix  form  as 
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— p» 


- 

n 

- 

^-2 

^  r,-i 

5"'' 

^0 

^-r,‘2 

^1 

0 

^  • 

■  ^  V3 

«2 

— 

0 

^P,  2 

■ 

0 

where  the  various  blocks  /?„,  R^,  and  /?,  are  given  by 


mO)  R{0,-\)  ••  R(0,-L^~P^^\) 

/?(0,1)  /?(0.0)  -  R(Q,-L^P^a2) 


R(OJL^^P^  l)  R{0,L^*P-2)  •  mO) 


Rik,-L,)  Rik,-L^-\)  •••  R{k,-L^-P^^l) 

R{k,-L^+l)  R{k,-L^)  ••  R{k,-L^-P^+2) 

/?(Jt,-Lj+P,-l)  R{k,-L^^P^~2)  ••  R{kfi) 

R{k,0)  R{k,-l)  ••  /?(*:,-/>,+ 1) 

^  R{k,l)  Rik,0)  ••  R{k,-P^*2)  (1.18) 

Rik,P^-2)  ..  /?(*:, 0) 

the  model  parameters  in  vector  form  are 

=  [fl(0,0)  am)  o{0,2)  ^7(04+P3-l)]^  (1-19) 

<7,  =  a(it^  +  l)  <7(^X2 +2)  -  fl(Uj+/’2-l)l’',  (1-2(11 

and 

5'”’  =  [a^  0  0  ••  0]'’. 

The  blocks  of  the  autocorrelation  matrix  have  three  structures  given  by  R^,  R,,,  and  R^ 

Except  for  the  upper  and  left  borders,  the  matrix  is  block-toeplitz  with  toeplitz  blocks 
Quarter-plane  support  is  a  special  case  of  NSHP  support  with  =  0. 


(1-17) 
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II.  SOLUTION  OF  NORMAL  EQUATION  USING  ITERATION 

The  solution  of  the  normal  equation  requires  inversion  of  the  autocorrelation 
matrix.  As  previously  stated,  direct  inversion  becomes  increasingly  computationally 
intensive  for  larger  matrices.  However,  the  toeplitz  block-toeplitz  structure  of  the 
autocorrelation  matrix  enables  the  inversion  of  the  matrix  via  an  iterative  method  which 
employs  successive  partitioning  of  the  normal  equation  [Ref.  81. 


A.  QUARTER-PLANE  SUPPORT 

To  develop  the  iterative  algorithm  we  first  divide  both  sides  of  (1-10)  by  This 
results  in  the  modified  normal  equation 


^0  ^-1  ^  2  ■■  ^  p,.l 

1 

1 _ 

S"”' 

/?,  /?„  /?_,  • 

0 

^1  ^0  ■■ 

Oi 

r: 

0 

/?p  1  Rp  _2  _3  ■■  ^0 

^r,-i 

0 

where  5'"’  -^[100-  0]^.  The  normal  equation  is  then  panitioned  as  follows. 
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(2-2) 


Equation  (2-2)  can  be  expressed  as 


where 

^0  ^-2  "  ^-/’,-2 

^1  ^0  ^-1  ” 

G,  =  /?,  /?„  •- 

^r,  2  ^p,-i  ^p,-*  "'  ^0 

<  -  ^P,2  ^.,3  •  ^,1- 

Y|  ”  t^O  ^1  ^2  ■■ 


(2-3) 


(2-4) 

(2-5) 

(2-6) 
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and 


(p,  =  11  0  0  ..  0]’". 


(2-7) 


Both  sides  of  (2-3)  are  then  premultiplied  by  a  matrix  F,  which  results  in 


where 


ff, 

Ro' 


(2-9) 


These  last  equations  suggest  an  iterative  solution.  In  particular,  (2-8)  may  be  expressed 
as 


p,-i  J 


(2-10) 


Equation  (2-10)  requires  the  inversion  of  the  submatrix  G,,  which  is  not  much  smaller 
than  the  autocorrelation  matrix,  and  hence  requires  nearly  as  many  multiplications  to 
invert.  Therefore,  G,,  y,,  and  tp,  may  be  further  partitioned  to  yield 
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(2-11) 


"2 

'  1 

_ 1 

<P2 

/// 

^0. 

^P,-2 

1 - 

0 

1 

where 


- 1 

>3 

0 

/?  , 

^-2  • 

/?, 

II 

^2 

^0  • 

^P,-3 

<n 

•  ^2 

(2-12) 


and 


H2  -  [Rp  .,  Rp.2  Rp,.,  ••  ^zl’ 

(2-13) 

r  T  T  T  T 

Y2  =  K  ^2  - 

(2-14) 

(p,  =  [1  0  0  ••  ojr 

(2-15) 

Premultiplication  of  both  sides  of  (2-11)  by  the  matrix  results  in 


Y2 -^2 '[92-^2 

Qp  2  ~  ~^0  ^2  y  2 


(2-16) 
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where 


w"  Ro' 

An  iterative  solution  to  (2-16)  can  then  be  formulated  as 


(2-17) 


(2-18) 


Since  (2-18)  requires  the  inversion  of  submatrix  Gj,  which  may  be  large,  the  matrices  G^, 
Yj,  and  <pj  are  further  partitioned.  Repetitive  partitioning  of  the  normal  equation  finally 
results  in 


where  Gp  .,  =  R„,  =  7?,,  Y^  .,  =  <pp  ,  =  5"”.  Premultiplication  of 

(2-19)  by  the  matrix  Fp  ,  yields  an  iterative  solution 


=  af-"] 


(2-20) 
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where 


(2-21) 


Combining  all  of  the  iterative  solutions  from  the  successive  levels  of  partitioning  results 
in  the  final  algorithm 


p.-t 


(1-1) 


<^l 


(1-1) 


l^n 


(2-22) 


where  j,  are  PjXl  vectors  of  AR  parameters,  af'^  =  R^'R.a^^  j  = 

1,2,. and  k  is  the  index  of  iteration.  Solution  of  (2-22)  requires  0(P]pIK) 

multiplications  where  K  is  the  number  of  iterations  required  for  convergence  to  the  true 
parameters. 

B.  NONSYMMETRIC  HALF-PLANE  SUPPORT 

As  noted  previously,  the  autocorrelation  matrix  which  results  from  NSHP  support 
consists  of  three  different  types  of  blocks.  The  derivation  of  the  iteration  follows  the 
same  procedure  of  succesive  partitioning  of  the  nonnal  equation  as  that  for  the  QP 
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support.  However,  the  assymmetry  of  the  autocorrelation  matrix  results  in  a  slightly 
different  final  recursion  given  by 


p.-i 


do  -  do  Aq 


(*-i) 


(t-i) 


/*! 

>*/ 


(2-23) 


where  a,  are  P^xl  vectors  of  AR  parameters,  Jo”’  =  ^o''^'”’.  J=-  1-2,..., 

Z*,-!,  and  k  is  the  index  of  iteration.  As  with  (2-22),  the  solution  of  (2-23)  requires 
0{P]pIK)  mutiplications  to  converge  to  the  true  values  of  the  parameters. 
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in.  COVARIANCE  METHOD 


The  iterative  solutions  to  the  normal  equations  given  in  (2-22)  and  (2-23)  assume 
that  the  autocorrelation  matrix  is  known  and  that  it  is  toeplitz  block-toeplitz.  In  general, 
however,  the  autocorrelation  matrix  must  be  estimated  from  the  available  data.  Assuming 
a  2-D  AR  model,  the  covariance  method  provides  a  means  to  estimate  the  autocorrelation! 
matrix  and  may  be  formulated  in  terms  of  linear  prediction. 

In  the  2-D  linear  prediction  problem,  the  error  between  the  true  value  of  a  random 
process  3nd  the  estimated  value  is  given  by 


(3-1) 

<’(«,, Hj)  =  •v(/J,,/J2)+ 53  53  . 

(1.;!  €  A 

(i.j)  *  (0.0) 

(3-2) 

The  objective  of  the  covariance  method  is  to  minimize  the  sum  of 

the  actual  squared 

errors  from  a  particular  set  of  these  terms  [Ref.  1].  Let  a(0,0)  =  1. 

written  as 

Then  (3-2)  can  be 

(i  .;■)  e  A 

This  last  equation  can  be  represented  in  vector  form  as 

(3-3) 

e  -  Xa 

(3-4) 
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where  the  rows  of  the  matrix  X  are  the  points  of  the  random  process  within  the  region  A 
as  the  filter  is  moved  over  the  data.  In  the  covariance  method  the  filter  is  positioned  so 
that  it  never  falls  even  partially  outside  of  the  available  data.  The  elements  of  the  vector 
e  are  the  error  between  the  observed  and  the  predicted  values  of  the  random  process  for 
each  position  of  the  filter.  Premultiplication  of  both  sides  of  (3-4)  by  and  application 
of  the  the  orthogonality  principle  results  in  the  following  normal  equation 

Ra  -  5,  (3-5) 

where  R  =  X^X,  S  =  X'^e  =  [5""^  0  ••  OJ^.and  S'®*  =  [o^  0  -  0]^.  The  structure  of 

R  depends  upon  the  shape  of  the  region  A.  Quarter-plane  and  NSHP  support  are 
discussed  below. 


A.  QUARTER-PLANE  SUPPORT 

If  the  region  A  is  a  quarter-plane  in  the  first  quadrant,  the  matrix  X  becomes 


X  - 


X 


v  • 


(3-6) 
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where  each  row  is  given  by 


x{n^^p,n^+q) 

x{n^*p-l,n^^q) 

x{ii,*p-\,n^+q-P^xl) 

.v(;/,  ip-P,  •  l,n^iq) 
.r(«,+/?-P,  I  \  ,n^*q-P^i  1) 


(3-7) 


and  N,  and  are  the  niiinber  of  columns  and  rows  in  the  2-D  random  process, 
respectively,  and  P,,  P^  are  the  dimensions  of  the  region  A  as  shown  in  Figure  4. 


Figure  4.  QP  support  used  to  estimate 
autocorrelation  matrix 
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The  sample  autocorrelation  matrix  then  becomes 


R  = 


^0.0 

^0.1 

^1,0 

(3-8) 


where  ^  matrix  R  is  symmetric;  however,  it  is  not  toeplitz  block-tocpiitz. 

Nevertheless,  an  iterative  solution  to  estimate  the  AR  parameters  may  still  be  obtained 
as  follows  [Ref.  8]. 

The  first  step  is  to  average  the  blocks  along  the  main  diagonal  of  the  autocorrelation 
matrix.  Refer  to  this  block  as  From  a  toeplitz  approximation  T  of  the  diagonal 
blocks  is  formed.  One  approximation  is  given  in  Ref.  9,  where  the  diagonal 
elements  of  T,  t(i).  are  given  by 


where  i  =  0,1....F2'1- 


'('» =  j—  E 


(3-9) 


2  J- 


The  normal  equation  is  then  successively  partitioned  as  shown 


previously. 

The  first  partitioning  of  the  normal  equation,  suniliar  to  that  in  (2-3),  leads  to 


(4- !  U 


.(*) 


^p,-\  -  Y 1 


(3-10) 


20 


The  block  may  be  expressed  as 


where  D^,  ,  is  the  difference  between  the  toeplitz  approximation  T  and 

Substituting  (3-11)  into  (3-10)  leads  to 


(3-11) 

p 


a^p[,^-T-'D 


(3-12) 


Subsequent  partitioning  of  the  normal  equation  then  results  in 


Yp,-i  ^p^-\  1 

^1  ~  ~"iA^F,-ty  p,-\ 


(*- 1  )1 


(3-13) 


where  y p  =  a^,Gp  ^  =  and  =  R^^.  The  submatrix  /?, ,  can  be  expressed 


as 


(3-14) 


and  /?„„  can  be  expressed  as 


^0.0  =  ^  +  ^0.0  • 


(3-15) 


Substituting  (3-14)  and  (3-15)  into  (3-13)  gives 
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(3-16) 


ar'] 

a<*>  =  -r-'D,  . 


Finally,  combining  all  of  the  succesive  iterations  leads  to 


a„  =ao  -T  D^^Oo  -T  ^  R^_,  a, 


aj=-T-'D  /?  V 

J  j.j  >  ' 


(3-17) 


where  a^^''  -  r  af^  -  and  j  =  1,2,... Equation  (3-17)  is  then 

iterated  to  solve  for  the  model  parameters  in  (3-5). 

B.  NSHP  SUPPORT 

An  iteration  similiar  to  that  used  for  QP  support  may  be  used  for  NSHP  support. 
Use  of  NSHP  suppoa  to  estimate  the  autocorrelation  matrix  results  in  a  matrix  with  an 
asymmetric  structure  similiar  to  that  of  (1-15).  The  sample  autocorrelation  matrix  can  be 
expressed  as 


i> 

0.0 

^U.l 

^or,- 1 

P 

'*.0 

p 

(3-18) 
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The  dimensions  of  the  blocks  along  the  top  and  left  edges  are  reduced  by  a  number  of 
rows  and  columns,  respectively,  equal  to  the  number  -L^  associated  with  region  A. 

The  development  of  the  iteration  begins  by  forming  a  toeplitz  approximation  T  of 
the  average  of  the  blocks  along  the  main  diagonal  of  R.  The  upper-left  block  is  not 
included  in  the  average,  however,  because  it  is  not  of  the  same  dimension  as  the  other 

blocks  along  the  diagonal.  A  separate  toeplitz  approximation  f  is  determined  for  this 
block. 

The  nonnal  equation  is  succesively  partitioned  as  discussed  previously  in  sections 
II. A  and  II. B  At  each  stage  of  the  partitioning,  the  blocks  along  the  main  diagonal  are 
expressed  as  the  sum  of  the  toeplitz  approximation  and  a  difference  matrix.  The  final 
fonn  of  the  algorithm  is 


.(*)  .(01  rf.-\  fZ  .(*-1)  T-IV'  O 

-T  -T  2^R„,a, 


(^I  T- - 1  ri  T*- I  D  'T'-l  %  ^  wy 

a,  =  -T  D  a,  -T  R  ..Or,  -T  '  }  R  a, 

)  J.)  )  J.o  ”  ^  ' 


(3-19) 


where  J,'”’  =  f  '5'"', 


'/?  X  ,  andy=  1,2,...P,-1. 


C.  COMBINED-QUADRANT  METHOD 

In  deriving  the  iteration  for  the  quarter-plane  case,  support  in  the  first  quadrant  was 
assumed.  However,  support  in  any  of  the  quadrants  may  have  been  used  without  any 
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loss  in  generality.  Due  to  the  hermitian  symmetry  and  toeplitz  block-toeplitz  property  of 
the  autocorrelation  matrix,  support  in  the  third  quadrant  results  in  an  identical  spectral 
estimate  to  that  obtained  from  support  in  the  first  quadrant.  In  general,  however,  spectral 
estimates  derived  from  support  in  the  first  and  second  quadrants  are  different. 

Spectral  estimates  obtained  from  quarter-plane  support  often  result  in  contours  of 
constant  pov/er  spectral  density  that  are  elliptical.  This  means  that  the  frequency  estimate 
in  one  direction  is  more  accurate  than  the  estimate  for  the  other  direction.  To  mitigate 
this  problem,  the  following  combined-quadrant  (CQ)  spectral  estimate  has  been  suggested 
[Ref.  10] 


where  P ^  and  P ^  are  the  spectral  estimates  obtained  from  QP  support  in  the  first  and 


second  quadrants,  respectively. 

It  can  easily  be  shown  that  QP  support  in  the  second  quadrant  results  in  the 
following  nomial  equation 


(3-21) 
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where  bj  =  [b{k,0)  bikA)  b{k,2)  ^(0,0)  =  1.  and 

5'“*  =  [a^  0  -  O]^  .  The  autocorrelation  matrix  is  identical  to  that  obtained  using  QP 

support  in  the  first  quadrant.  However,  the  AR  parameters  fc,  are  not  the  same  in  general. 

The  estimation  of  the  AR  parameters  is  obtained  by  an  iteration  similiar  to  (3-17). 
The  difference  is  the  manner  in  which  the  normal  equation  is  partitioned.  Rather  than 
beginning  the  partitioning  at  the  lower-right  comer  and  continuing  toward  the  opposite 
comer,  as  previously  shown  with  support  in  the  first  quadrant,  the  partitioning  begins  with 
the  upper-left  comer  and  continues  to  the  opposite  comer.  The  toeplitz  approximation  T 
is  the  same  as  that  used  for  support  in  the  first  quadrant.  The  iteration  is  therefore 


1^1 


(3-22) 


where  -  T '5'"’,  bf' 


r  'Rp  ,  j  =  1,  2, ....  P,-l.  The  CQ  method  requires 


2  2 

0(2P,P2K)  multiplications  to  converge  to  the  true  parameter  values. 


IV.  RESULTS  OF  ITERATIVE  SPECTRAL  ESTIMATION 


The  performance  of  the  iterative  spectral  etimation  method  was  investigated  by 
testing  its  ability  to  detect  multiple  sinusoids  in  white  noise.  The  data  was  given  by 

x{n^,n2)  - B jCos{Q.l25ii^+0.l25n^)+B^cosi0.33n^+033n^)+w{n^,n2) 

where  are  the  amplitudes  of  the  sinusoids,  and  wini.n^)  is  a  sample  function  of  zero 
mean  white  noise  with  unit  variance.  The  signal-to-noise  ratio  (SNR)  of  the  random 
process  is  defined  by  (Ref.  1] 

SNR  (4-2) 


where  m  is  the  number  of  sinusoids  present. 

Two  sizes  of  data  sets  were  used  to  evaluate  the  performance  of  the  iterative 
spectral  estimation  method:  16x16  and  8x8.  For  each  size  a  SNR  of  10  dB  and  0  dB 

were  used.  The  SNR  was  altered  by  varying  the  values  of  B,  while  holding  G^  constant. 

Comparisons  were  made  between  QP  and  NSHP  support  and  the  CQ  method. 

Spectra]  estimates  of  the  16x16  data  set  with  a  10  dB  SNR  are  shown  in  Figure  5. 
Best  results  were  obtained  for  QP  support  after  12  iterations  for  P,  =  Fj  =  4.  The 
estimated  frequencies  do  not  match  the  true  values,  and  there  are  many  spurious  peaks. 
Other  values  of  P,  and  Pj  faded  to  produce  more  accurate  results. 
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Figure  5.  Spectral  estimates  of  16x16  data 


Better  results  were  obtained  using  NSHP  support  with  f’,=P2=4,  and  -  2  .  After 

one  iteration  an  acceptable  spectral  estimate  was  obtained.  The  estimated  frequencies  are 
closer  to  the  true  values,  and  the  spectral  peaks  are  sharp. 

The  CQ  method  produced  even  better  results.  An  acceptable  spectral  estimate  was 
obtained  after  a  single  iteration  with  P,  =  Pj  =  3.  The  spectral  peaks  are  sharp,  and  the 
estimated  frequencies  are  closer  to  the  true  values  than  with  NSHP  support. 

Results  of  spectral  estimation  applied  to  a  0  dB  SNR  data  set  of  size  16x16  are 
shown  in  Figure  6.  Meaningful  spectral  estimates  could  not  be  obtained  using  QP  support 
for  any  values  of  P,  and  Pj.  Employment  of  NSHP  support  led  to  a  useful  spectral 
estimate  after  three  iterations  for  P,=P2=4,  and  £2= -2.  However,  the  spectral  peaks  are 
broad. 

The  CQ  method  yielded  a  useful  spectral  estimate  after  one  iteration  for  P,  and  P2 
equal  to  three.  The  estimated  frequencies  are  very  close  to  the  true  values,  but  the 
spectral  peaks  are  very  broad. 

Spectral  estimates  of  a  10  dB  SNR  data  set  of  size  8x8  are  seen  in  Figure  7.  After 
five  iterations  QP  support  produced  a  mediocre  spectral  estimate  for  both  P,  and  P2  equal 
to  four.  The  estimated  frequencies  are  very  close  to  the  true  values.  However,  the 
presence  of  spurious  peaks  makes  the  spectral  estimate  very  difficult  to  interpret 
accurately,  in  spite  of  the  sharp  spectral  peaks. 

A  more  nearly-accurate  spectral  estimate  was  obtained  from  NSHP  support  after  a 
single  iteration  for  P^=Pj=4,  and  L2=  -  2.  The  accuracy  of  the  higher-frequency  estimate 
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8x8  data 


is  very  good,  whereas  the  resolution  of  the  lower-frequency  estimate  is  only  fair. 

The  CQ  method  produced  a  very  accurate  spectral  estimate  after  one  iteration  for 
/’,  and  p2  equal  to  three.  As  seen  in  the  figure,  the  spectral  peaks  are  very  sharp,  and  the 
estimated  frequencies  are  very  close  to  the  true  values. 

Only  the  CQ  method  produced  an  acceptable  spectral  estimate  of  the  0  dB  SNR  data 
set  of  size  8x8.  A  meaningful  spectral  estimate  was  obtained  after  a  single  iteration  for 
R,  and  equal  to  three,  as  seen  in  Figure  8.  The  frequency  estimates  are  close  to  the  true 
values,  but  the  spectral  peaks  are  broad. 

Several  other  cases  were  tested  for  different  smusoidal  frequencies  in  white  noise. 
Additionally,  the  values  for  the  SNR  were  varied.  In  all  of  these  cases,  the  results  were 
consistent  with  those  presented  above. 
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V.  CONCLUSIONS 


The  results  presented  in  Chapter  FV  show  that  high-resolution  spectral  estimates  can 
be  obtained  using  iterative  methods.  In  several  cases,  a  single  iteration  was  sufficient  to 
produce  meaningful  results.  The  net  gain  is  a  reduction  in  the  number  of  computations 
required  to  estimate  the  spectra. 

Comparisons  of  spectral  estimates  obtained  from  QP  and  NSHP  support  and  the  CQ 
method  show  that  the  latter  produced  superior  results.  Although  the  CQ  method  requires 

OilP'lPl)  multiplications  to  arrive  at  a  solution  as  compared  to  0{P^^Pl)  for  QP  and 

NSHP  support,  it  requires  a  smaller  region  of  support.  Tlierefore,  the  CQ  method  not  only 
provides  spectral  estimates  of  superior  resolution  to  QP  and  NSHP  support,  it  also 
accomplishes  these  results  in  fewer  iterations. 

It  was  observed  that  the  resolution  of  the  spectral  estimates  did  not  necessarily 
improve  with  more  iterations.  In  the  case  of  the  CQ  method,  the  resolution  did  not 
unprove  significeinlly  after  the  first  iteration.  Results  for  QP  and  NSHP  support  indicate 
that  the  estimates  of  the  AR  parameters  converge  to  the  true  values,  and  may  then 
diverge.  This  phenomenom  is  not  understood  and  provides  a  basis  for  continued  study. 

Additional  areas  of  research  arise  from  the  findings  of  this  thesis.  Tire  fast  solution 
of  the  nonnal  equation  via  iteration  leads  to  the  possibility  of  following  spectral  lines  of 
slowly  time-varying  signals.  Because  the  values  of  the  AR  parameters  are  dependent  upon 
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previous  values,  the  previous  values  could  possibly  be  used  to  re-initialize  the  estimates 
when  the  random  process  changes  in  time.  Exploration  of  this  idea  would  serve  as  good 
follow-on  study.  Another  follow-on  study  is  the  posible  application  of  iterative  methods 
to  the  multichannel  case,  where  the  AR  paraneters  are  in  matrix  form. 
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